[file path] = uigetfile('*metadata.mat');

load([path file]);

cd(path);

vidList = vidList(logical([vidList.valid])); %Remove any non-valid videos
%%
vis = 2;  %Viusalize to console?
mov = 1; %Make output videos?

trialMetadata = struct();

for i = 1:numel(vidList)
    %% Load video and calculate background
    
    vid = VideoReader(vidList(i).name);
    
    bg = calculateBackgroundFromVideo(vid, 20, 0.8); %This step takes an incredibly long time bc of indexing the huge videos
    %% Open output video for visualization?
    
    if mov == 1
        outvid = VideoWriter(strrep(vid.name, '.avi', '_sampleTracking.avi'));
        open(outvid);
    end
    
    %% Get video metadata
    leftBaitCentroid = regionprops(vidList(i).leftBaitPoly);
    leftBaitCentroid = leftBaitCentroid.Centroid;
    
    rightBaitCentroid = regionprops(vidList(i).rightBaitPoly);
    rightBaitCentroid = rightBaitCentroid.Centroid;
    
    lbPoly  = vidList(i).leftBaitPoly;
    rbPoly = vidList(i).rightBaitPoly;
    
    lbPolyPoints = mask2poly(lbPoly);
    rbPolyPoints = mask2poly(rbPoly);
    %% Set up cross-frame loop
    nframes = vid.NumberOfFrames;
    
    nsamples = 4000;
    
    permissableSizeRange = [150 800]; %Only look at blobs between 50 and 500 pixels in area
    samples = round(linspace(1,nframes, nsamples)); %Select a subset of frames
    
    %% Set up empty output variables
    outVars = {'frame', 'meanDistanceToLeftBait', 'meanDistanceToRightBait', 'PercentOnLeftBait', ...
        'PercentOnRightBait', 'meanDistanceToLeftBaitBox', 'meanDistanceToRightBaitBox'};
    outputData = nan(nsamples, numel(outVars));
    
    %%
    %for j = 1:nsamples
    for j = 1:50
        
        %%
        im = rgb2gray(read(vid, samples(j)));
        imd = double(im) - double(bg);
        
        %% calculate binary thresholded image
        thresh = -30;
        imb = imd < thresh;
        %imshow(imb);
        
        %Erode image
        se = strel('disk', 2);
        imb = imerode(imb, se);
        
        %Dilate image
        imb = imdilate(imb, se);
        
        %Calculate regionprops info on binary image
        rp = regionprops(imb, {'Area', 'Centroid', 'BoundingBox', 'Orientation', 'MinorAxisLength', 'MajorAxisLength'});
        rp = rp([rp.Area] > permissableSizeRange(1) & [rp.Area] < permissableSizeRange(2));
        arena = vidList(i).arenaPoly;
        
        cents = round(cat(1,rp.Centroid)); %Subset to
        for k = 1:numel(cents(:,1))
            cents(k,3) =  arena(cents(k,2), cents(k,1));
        end
        
        rp = rp(logical(cents(:,3)));
        
        
        
        
        %% calculate output metrics
        
        
        centroids = cat(1,rp.Centroid);
        
        %         imshow(im);
        %         hold on
        %         plot(centroids(:,1), centroids(:,2), 'ro');
        %         plot(leftBaitCentroid(1), leftBaitCentroid(2), 'bo');
        %         plot(rightBaitCentroid(1), rightBaitCentroid(2), 'bo');
        %         hold off
        
        
        
        
        %Distances
        
        distLeftBait = sqrt((centroids(:,1) - leftBaitCentroid(1)).^2 + (centroids(:,2) - leftBaitCentroid(2)).^2);
        distRightBait = sqrt((centroids(:,1) - rightBaitCentroid(1)).^2 + (centroids(:,2) - rightBaitCentroid(2)).^2);
        
        
        %Create empty data vectors
        onLeftBait = nan(numel(centroids(:,1)),1);
        onRightBait = nan(numel(centroids(:,1)),1);
        
        distToLeftBaitBox = nan(numel(centroids(:,1)), 1);
        distToRightBaitBox = nan(numel(centroids(:,1)), 1);
        
        for zz = 1:numel(onLeftBait) %Looop across centroids
            centCur = round(centroids(zz,:));
            
            %Distance to box
            distToLeftBaitBox(zz) = p_poly_dist(centCur(1), centCur(2), lbPolyPoints(:,1), lbPolyPoints(:,2));
            distToRightBaitBox(zz) = p_poly_dist(centCur(1), centCur(2), rbPolyPoints(:,1), rbPolyPoints(:,2));
            
            %Location on/off bait
            
            onLeftBait(zz) = lbPoly(centCur(2), centCur(1));
            onRightBait(zz) = rbPoly(centCur(2), centCur(1));
        end
        
        onLeftBait = logical(onLeftBait);
        onRightBait = logical(onRightBait);
        
        %Set points within bait box to zero
        distToLeftBaitBox(onLeftBait) = 0;
        distToRightBaitBox(onRightBait) = 0;
        
        
        
        
        
        
        if vis == 1
            imshow(im);
            hold on
            plot(centroids(:,1), centroids(:,2), 'r.', 'MarkerSize', 20);
            %plot(leftBaitCentroid(1), leftBaitCentroid(2), 'b.');
            %plot(rightBaitCentroid(1), rightBaitCentroid(2), 'b.');
            plot(centroids(onLeftBait,1), centroids(onLeftBait,2), 'go', 'MarkerSize', 15);
            plot(centroids(onRightBait,1), centroids(onRightBait,2), 'bo', 'MarkerSize', 15);
            xlim([0 1920]);
            ylim([0 1080]);
            
            hold off
            
        elseif vis == 2
            
            imshow(im);
            hold on
            plot(centroids(:,1), centroids(:,2), 'r.', 'MarkerSize', 20);
            plot(leftBaitCentroid(1), leftBaitCentroid(2), 'g.', 'MarkerSize', 20);
            plot(rightBaitCentroid(1), rightBaitCentroid(2), 'b.', 'MarkerSize', 20);
            
            for zz = 1:numel(centroids(:,1))
                plot([leftBaitCentroid(1) centroids(zz,1)], [leftBaitCentroid(2) centroids(zz,2)], 'g', 'LineWidth', 0.1);
                plot([rightBaitCentroid(1) centroids(zz,1)], [rightBaitCentroid(2) centroids(zz,2)], 'b', 'LineWidth', 0.1);
                
            end
            
            xlim([50 1200]);
            ylim([0 1000]);
            
            hold off
            % % Alternative - plot body ellipses
            %             imshow(im);
            %             alphamask(imb, [0 1 0], 0.5);
            %             t = linspace(0,2*pi,50);
            %
            %             hold on
            %             for k = 1:length(rp)
            %                 a = rp(k).MajorAxisLength/2;
            %                 b = rp(k).MinorAxisLength/2;
            %                 Xc = rp(k).Centroid(1);
            %                 Yc = rp(k).Centroid(2);
            %                 phi = deg2rad(-rp(k).Orientation);
            %                 x = Xc + a*cos(t)*cos(phi) - b*sin(t)*sin(phi);
            %                 y = Yc + a*cos(t)*sin(phi) + b*sin(t)*cos(phi);
            %                 plot(x,y,'r','Linewidth',1)
            %             end
            %             hold off
        end
        
        if mov == 1
            F = getframe(gcf);
            writeVideo(outvid, F);
        end
        %% Export analyzed data
        outputData(j,1) = samples(j);
        outputData(j,2) = mean(distLeftBait);
        outputData(j,3) = mean(distRightBait);
        outputData(j,4) = mean(onLeftBait);
        outputData(j,5) = mean(onRightBait);
        outputData(j,6) = mean(distToLeftBaitBox);
        outputData(j,7) = mean(distToRightBaitBox);
    end
    
    if mov == 1
        close(outvid);
    end
    
    %% export data to disk for this video
    outputData = array2table(outputData);
    outputData.Properties.VariableNames = outVars;
    vidName = vidList(i).name;
    vidFolder = vidList(i).folder;
    
    outfile = strrep(vidName, '.mov', '_tracked.csv');
    
    writetable(outputData, outfile);
    %%
    trialMetadata(i).vidName = vidName;
    trialMetadata(i).vidFolder = vidFolder;
    trialMetadata(i).trackedDatafile = outfile;
    trialMetadata(i).concentration = vidName(7);
    trialMetadata(i).leftBait = vidName(12);
    trialMetadata(i).rightBait = vidName(17);
    
    i
end

%% save trial metadata file

writetable(struct2table(trialMetadata), 'trialMetadata.csv');

%
% outstr = strrep(file, 'metadata.mat', 'metadata_analyzed.mat');
% save(outstr, 'vidList');